%This function unfolds the matrices reported by Dynare++ up to a fourth order approximation
%The assumption is that dynare++ .mat files are saved on the path, with '_#' indicating the 
%underlying approximation number.
%
% Our notation is as follows
% z_t = f([x_t-1;u_t])
% where 'x'    the state variables
%       'u'    the innovations 
% and v_t = [x_t-1;u_t]
%
% IMPORTANT: the codes are made for running DYNARE++ with the option --no-centralize
% By Martin M. Andreasen

function [f_ss,fv,fvv,fss,fvvv,fssv,fsss,fvvvv,fssvv,fsssv,fssss,sigma,nv,nu,nx,nz,dyn_vars,dyn_state_vars,dyn_shocks] ...
    = Dynarepp_Unfold_Matrices_4th(model,order_app);

%Load model to initialise counter and aux_vars
load([model,'_1','.mat']);

%Auxiliary definitions
nv      = size(dyn_state_vars,1);   % Number of state variables and innovations
nu      = size(dyn_vcov_exo,1);     % Number of innovations
nx      = nv-nu;                    % Number of state variables
nz      = size(dyn_vars,1);         % Number of control variables + state variables

% The cholesky factorisation of the covariance matrix
[sigma,error_mes] = chol(dyn_vcov_exo);       
if error_mes ~= 0
    sigma = diag(sqrt(diag(dyn_vcov_exo)));
    disp('Calculaions done for diagonal shock covariance matrix')
end

%Load matrices with coefficients:
%Deterministic steady state
f_ss = dyn_ss;
%f_ij contains the matrix dyn_g_j corresponding to an i-th order approximation
for i = 1:order_app
    load([model,'_',num2str(i),'.mat']);
    for j = 0:i
        eval(['f_',num2str(i),num2str(j),'=dyn_g_',num2str(j),';']);
    end
end

%The first order terms
fv             = f_11;

% The second order terms 
% We multiply by 2 because the terms in Dynare++ are already multiplied by 1/2 (see below)
if order_app > 1
    fss            = 2*f_20;
    % we unfold the matrix f_22
    fvv            = zeros(nz,nv,nv);
    index = 0;
    for alfa1=1:nv
        for alfa2=alfa1:nv
            index = index + 1;
            fvv(:,alfa1,alfa2) =  f_22(:,index);
            if alfa1 ~= alfa2
               fvv(:,alfa2,alfa1) =  f_22(:,index);
            end
        end
    end
    fvv = 2*fvv;
else
    fss = zeros(nz,1);
    fvv = zeros(nz,nv,nv);
end

% The third order terms are given by f_33 and by the difference between f_31 and f_11
% We multiply by 6 because the terms in Dynare++ are already multiplied by 1/6 (see below)
if order_app > 2
    fssv           = (6/3)*(f_31-f_11);            %the term f_31 is multiplied by 3/6 in Dynare ++ because there are
                                                   %3 combinations (fssz,fszs,fzss) in these third order terms. 
    % Fx. f_11 = gv and f_31 = gv+3/6*gssv. So, f_31-f_11 = 3/6*gssv, and therefore 6/3*(f_31-f_11) = gssv. 
    
    % We unfold the matrix f_33 and multiplt by 6
    fvvv           = zeros(nz,nv,nv,nv);
    index          = 0;
    index_terms    = 0;
    for alfa1=1:nv
        for alfa2=alfa1:nv
            for alfa3=alfa2:nv
                index = index + 1;
                fvvv(:,alfa1,alfa2,alfa3) = f_33(:,index);
                index_terms = index_terms + 1;
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %fvvv(:,alfa1,alfa1,alfa3)= fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa1,alfa3,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa3,alfa1,alfa1) = fvvv(:,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %fvvv(:,alfa1,alfa2,alfa2)= fvvv(:,alfa1,alfa2,alfa3);                
                    fvvv(:,alfa2,alfa1,alfa2) = fvvv(:,alfa1,alfa2,alfa3);                
                    fvvv(:,alfa2,alfa2,alfa1) = fvvv(:,alfa1,alfa2,alfa3);                                
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %fvvv(:,alfa1,alfa2,alfa3) = fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa1,alfa3,alfa2) = fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa3,alfa1,alfa2) = fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa3,alfa2,alfa1) = fvvv(:,alfa1,alfa2,alfa3);
                    fvvv(:,alfa2,alfa3,alfa1) = fvvv(:,alfa1,alfa2,alfa3); 
                    fvvv(:,alfa2,alfa1,alfa3) = fvvv(:,alfa1,alfa2,alfa3);                 
                end
            end
        end
    end
    fvvv = 6*fvvv;
    fsss = zeros(nz,1);    %Dynare++ solves only for Gaussian shocks
else
    fsss = zeros(nz,1);
    fssv = zeros(nz,nv);
    fvvv = zeros(nz,nv,nv,nv);
end


% The fourth order terms are given by f_44, f_40, the difference between f_42 and f_22
% We multiply by 6 because the terms in Dynare++ are already multiplied by 1/24 (see below)
if order_app > 3
    fssss           = (24/1)*(f_40-f_20);            
    % Fx. f_20 = 1/2*gss and f_40 = 1/2*gss+1/24*gssss. So f_40 - f_20 = 1/24*gssvv, and therefore gssss = 24/1*(f_40 - f_20) 
    fsssv           = zeros(nz,nv);                %Dynare++ solves only for Gaussian shocks
    fssvv_tmp       = (24/6)*(f_42-f_22);          %the term f_42 is multiplied by 6/24 in Dynare ++ because there are
                                                   %6 combinations in these fourth order terms

    % Fx. f_22 = 1/2gvv and f_42 = 1/2gvv+6/24gssvv. So, f_42-f_22 = 6/24gssvv, and therefore 24/6*(f_42-f_22) = gssvv.                                                       
    fssvv           = zeros(nz,nv,nv);
    index = 0;
    for alfa1=1:nv
        for alfa2=alfa1:nv
            index = index + 1;
            fssvv(:,alfa1,alfa2) =  fssvv_tmp(:,index);
            if alfa1 ~= alfa2
               fssvv(:,alfa2,alfa1) = fssvv_tmp(:,index);
            end
        end
    end
    
    % We unfold the matrix f_44 and multiply by 24
    fvvvv           = zeros(nz,nv,nv,nv,nv);
    %fvvvv           = NaN(nz,nv,nv,nv,nv);    
    index           = 0;
    index_terms     = 0;
    for alfa1=1:nv
        for alfa2=alfa1:nv
            for alfa3=alfa2:nv
                for alfa4=alfa3:nv
                    index = index + 1;
                    fvvvv(:,alfa1,alfa2,alfa3,alfa4) = f_44(:,index);
                    index_terms = index_terms + 1;

                    %1) Three identical and one different, for instance(1,1,1,2)
                    if alfa1 == alfa2 && alfa1 == alfa3 && alfa1~= alfa4
                        index_terms = index_terms + 3;
                        %fvvvv(:,alfa1,alfa1,alfa1,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa1,alfa4,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa4,alfa1,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                        
                        fvvvv(:,alfa4,alfa1,alfa1,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                    end
                    %2) Three identical and one different, for instance(1,2,2,2)
                    if alfa1 ~= alfa2 && alfa2 == alfa3 && alfa3 == alfa4
                        index_terms = index_terms + 3;
                        %fvvvv(:,alfa1,alfa2,alfa2,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa1,alfa2,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa2,alfa1,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa2,alfa2,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                        
                    end               
                    %3) Two identical and two different, for instance(1,1,2,2)
                    if alfa1 == alfa2 && alfa3 == alfa4 && alfa1 ~= alfa3
                        index_terms = index_terms + 5;
                        %fvvvv(:,alfa1,alfa1,alfa3,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa3,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa3,alfa1,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                        
                        fvvvv(:,alfa1,alfa3,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);   
                        fvvvv(:,alfa3,alfa1,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                        
                        fvvvv(:,alfa3,alfa3,alfa1,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                        
                    end
                    %4) To identical but three different, for instance (1,2,3,3) 
                    if alfa1 ~= alfa2 && alfa3 == alfa4 && alfa1~= alfa3 && alfa2~=alfa3
                        index_terms = index_terms + 11;
                        %fvvvv(:,alfa1,alfa2,alfa3,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa1,alfa3,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa3,alfa2,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa3,alfa3,alfa1,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa1,alfa3,alfa3,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa3,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa2,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa3,alfa1,alfa3,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa2,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa3,alfa1,alfa2,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa1,alfa3,alfa2,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa3,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);                      
                    end
                    %5) Two identical but three different, for instance (1,1,2,3)
                    if alfa1 == alfa2 && alfa3~= alfa4 && alfa1~=alfa3 && alfa2 ~= alfa3
                        index_terms = index_terms + 11;  
                        %fvvvv(:,alfa1,alfa1,alfa3,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa1,alfa4,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa4,alfa1,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa4,alfa3,alfa1,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa1,alfa4,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa4,alfa1,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa1,alfa3,alfa4,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa4,alfa3,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa1,alfa3,alfa1,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa4,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa3,alfa1,alfa1,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa4,alfa1,alfa1,alfa3) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                    end
                    %6) Two identical but three different, for instance (1,2,2,3)
                    if alfa2 == alfa3 && alfa1 ~= alfa4 && alfa1 ~= alfa2 && alfa3 ~= alfa4
                        index_terms = index_terms + 11;  
                        %fvvvv(:,alfa1,alfa2,alfa2,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa4,alfa2,alfa2,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa2,alfa4,alfa2,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa1,alfa2,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa2,alfa4,alfa1,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa1,alfa4,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa4,alfa2,alfa1,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa2,alfa4,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa4,alfa1,alfa2,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa1,alfa4,alfa2,alfa2) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        
                        fvvvv(:,alfa2,alfa2,alfa1,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        fvvvv(:,alfa2,alfa2,alfa4,alfa1) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                    end
                    %7) all different
                    if alfa1 ~= alfa2 && alfa1 ~= alfa3 && alfa1 ~= alfa4 && alfa2 ~= alfa3 && alfa2 ~= alfa4 && alfa3 ~= alfa4
                        index_terms = index_terms + 24;
                        %fvvvv(:,alfa1,alfa2,alfa3,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        % Note this implementation implies that we
                        % "reproduce fvvvv(:,alfa1,alfa2,alfa3,alfa4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4)"
                        % This implies that index_terms > nz^4
                        idx = perms([alfa1,alfa2,alfa3,alfa4]);
                        for i=1:24
                            k1 = idx(i,1);
                            k2 = idx(i,2);
                            k3 = idx(i,3);
                            k4 = idx(i,4);
                            fvvvv(:,k1,k2,k3,k4) = fvvvv(:,alfa1,alfa2,alfa3,alfa4);
                        end
                    end 
                end
            end
        end
    end
    % Check if any of the terms in fvvvv are not assigned a real number (for debugging)
    %any(any(any(any(any(isnan(fvvvv))))))
    fvvvv = 24*fvvvv;
else
    fssss = zeros(nz,1);
    fsssv = zeros(nz,nv);
    fssvv = zeros(nz,nv,nv);    
    fvvvv = zeros(nz,nv,nv,nv,nv);
end

% Change format for dyn_vars, dyn_state_vars, dyn_shocks so they are cells
dyn_vars = cellstr(dyn_vars);
dyn_state_vars = cellstr(dyn_state_vars);
dyn_shocks = cellstr(dyn_shocks);
end
 


% ************************************************************************
% For the sake of completness, we here reproduce the way the information is
% stored in Dynare++ (a reply by Ondra K.)
%
% The rows of the matrices are clear, right? The ordering of rows is given by dyn_vars vector. 
% Columns determine with what deviations the number in the column is to be multiplied. 
% NOTE THAT the deviations are not from steady state, the deviations are from a fix-point, 
% which is saved as dyn_ss. So this means that you cannot find dyn_g_0, since it is zero, 
% since dyn_ss is chosen so that dyn_g_0 is zero (from the definition of the fix-point). 
% 
% The matrices are already mulitplied with the multipliers 1/2, 1/6 etc. 
% 
% The columns are ordered as folded indices of variables stored in the vector dyn_state. 
% You can find more details first in kord.pdf and then in tl.pdf downloadable from the documentation page. 
% Everything is described there. However, I provide a small example. Let a model have the states [K,C,EPS]. 
% If you assign indices to these variables in increasing order, say [K=0,C=1,EPS=2], 
% then the folded indices will be all growing sequences consisting from 0,1,2 in alphabetical ordering. 
% This is, for the third order: 
% 000 
% 001 
% 002 
% 011 
% 012 
% 022 
% 111 
% 112 
% 122 
% 222 
%
% for the fourth order: 
% 0000 
% 0001 
% 0002 
% 0011 
% 0012 
% 0022 
% 0111 
% 0112 
% 0122 
% 0222 
% 1111 
% 1112 
% 1122 
% 1222 
% 2222 
%
